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Strong coupling of finite element methods for the 

Stokes-Darcy problem* 
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Abstract 



The aim of this paper is to propose a systematic way to obtain convergent finite 

P^ element schemes for the Darcy-Stokes flow problem by combining well-known mixed 

finite elements that are separately convergent for Darcy and Stokes problems. In the 

^ approach in which the Darcy problem is set in its natural H(div) formulation and 



the Stokes problem is expressed in velocity-pressure form, the transmission condition 
ensuring global mass conservation becomes essential. As opposed to the strategy 
that handles weakly this transmission condition through a Lagrange multiplier, we 
impose here this restriction exactly in the space of global velocity field. Our analysis 
of the Galerkin discretization of the resulting problem reveals that, if the mixed finite 
element space used in the Darcy domain admits an H(div)-stable discrete lifting of 
the normal trace, then it can be combined with any stable Stokes mixed finite 
element of the same order to deliver a stable global method with quasi-optimal 
convergence rate. Finally, we present a series of numerical tests confirming our 



r^ theoretical convergence estimates. 

cn 

(Ni 1 Introduction 

>■ In this paper we are interested by the mixed finite element approximation of the coupled 

k> Darcy-Stokes problem. The Darcy-Stokes coupled system provides a linear model for the 

J-j simulation of incompressible flows in heterogeneous media. It is governed by the Stokes 

equations in one part of the domain, while in the other part, the flow is described by 
a standard second order elliptic equation derived from Darcy's law and conservation of 
mass. Proper transmission conditions must also be prescribed on the boundary common 
to the two media: conservation of mass enforces continuity of the normal velocities at this 
interface and conservation of momentum enforces the balance of the normal stresses. A 
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further interface condition, supported by empirical evidence and known as the Beavers- 
Joseph- Saff man condition, must also be taken into account, cf. [2l |20| [T6]. 

The development of suitable numerical methods for the Darcy-Stokes flow interaction 
has become a subject of increasing interest during the last decade. Discacciati et al. [7] 
provided the first theoretical study of the problem. The Galerkin scheme discussed in [7j 
is based on a standard finite element method for second order elliptic problems in the 
Darcy domain. Motivated by the need of obtaining direct finite element approximations 
of the Darcy flux, most of the approaches (cf. for instance [HI [3l [TOl [TTl [12] ) consider 
now an H(div)-flux conforming formulation in the Darcy domain. Hence, in this case, 
a mixed variational formulation in the porous media is coupled with the usual velocity- 
pressure formulation in the Stokes domain. We also point out that other finite element 
discretization strategies have been considered for this problem. For instance, a discon- 
tinuous Galerkin method is described in [19] and a stabilized finite element method is 
presented in [6]. 

In the mixed approach, the transmission condition that guaranties the equilibrium of 
the normal velocities becomes essential. One approach to deal with this constraint consists 
in enforcing it weakly by means of a Lagrange multiplier representing the trace of the 
Darcy pressure on the common interface, see for example [TSl [101 112] • One known feature 
of the finite element discretization of this formulation is that two meshes are required on 
the transmission boundary satisfying a stability condition between their corresponding 
mesh sizes, which certainly constitutes a quite cumbersome restriction. 

In this paper we are interested in the strong coupling that considers the inclusion of 
the essential transmission condition directly in the definition of the space to which the 



Darcy flux and the fluid velocity belong (see (2.9) below). This formulation prevents 
from introducing a further unknown (the Lagrange multiplier) simplifying by the way the 
saddle point structure of the problem. Two finite element discretizations have already 
been proposed for this formulation of the problem. In the first one [17] the authors 
provide a unified approximation in the whole domain giving rise to a global H(div)- 
conforming scheme that is nonconforming in the fluid domain. In the second one [3], a 
mortar discretization that combines the Bernardi-Raugel element and the lowest order 
Raviart-Thomas element is studied. 

Our concern in this paper is to figure out how to combine any of the common mixed 
finite element methods for Darcy and Stokes problems to end up with a stable scheme for 
the coupled problem. We consider stable conforming Galerkin schemes in each subdomain 
and establish the rules under which these schemes can be matched to form a global 
convergent method. We will see that the only requirement for the viability of such a 
coupling is the existence of a discrete stable lifting of the normal trace in H(div). We can 
prove that such a lifting is available in the bi-dimensional case for the Raviart-Thomas and 
the Brezzi-Douglas-Marini elements while, in the three-dimensional case, to obtain the 
same result, we need to assume that the triangulation is quasi-uniform in a neighborhood 
of the transmission boundary. 

The paper is organized as follows. In Section |2] we decribe the governing equations, 
derive a mixed variational formulation of the problem and show that it is well-posed. Our 
Galerkin scheme is introduced in Section |3] by means of abstract finite elements spaces in 
each subdomain. We also establish in this section the conditions guarantying the stability 



of our approximation method. In Section [4], we discuss the approximation properties of 
the discrete space for the global velocity field and provide our main convergence result. 
The key hypotheses providing the convergence of our numerical scheme are discussed in 
Sections |5] and [6j We show that they are satisfied for the most common mixed finite ele- 
ments for Darcy and Stokes problems under mild conditions on the family of finite element 
triangulations. In Section [7] we provide some examples obtained by combining well-known 
finite elements for the Stokes problem, such as the MINI element and the Bernardi-Raugel 
element, with the Raviart-Thomas and the Brezzi-Douglas-Marini elements in order to 
obtain a global scheme for the coupled Darcy-Stokes problem. Finally, in Section [8} we 
present a set of numerical experiments that confirm the converge rate predicted by the 
theory of the examples described in the former section. 

Notation and background. Boldface fonts will be used to denote vectors and vector 
valued functions. Also, if if is a vector space of scalar functions, H will denote the space 
of M*^ valued functions whose components are in H, endowed with the product norm. 

Given an integer m > 1 and a bounded Lipschitz domain O C M'', {d = 2,3), we 
denote by || • ||_H-m(c)) the norm in the usual Sobolev space H"^{0), cf. p. For economy 
of notation, (-, ■)o stands for the inner product in L^(0) and || ■ Ho is the corresponding 
norm. We recall that H^^'^(dO) represents the image of H^{0) by the trace operator. 
Its dual with respect to the pivot space L'^{dO) is denoted H~^/'^{dO). If S is a part of 
the Lipschitz boundary dO^ we denote by i^gg (S) the space of functions from H^f^ijl) 
whose extension by zero to the whole dO belongs to H^^'^{dO). We will consider the dual 
i^QQ (E) of i^QQ (E) with respect to the pivot space -^^(E) and the angled bracket (■, ■)^ 
will be used for the -^^(E) inner product and its extension as the duality product between 
Hqq (E) and H^q (E). For definition and basic properties of the space H(div, (9), we 
refer to [11]. We will denote by Ho(div, (9) the subspace of fields from H(div, (9) with 
zero normal trace on the boundary dO. 

Given an open set O or the closure of an open set, Pa:(0) will denote the space of 
(i— variate polynomials of degree not greater than k. The same notation will be applied for 
polynomials defined on flat {d — l)-dimensional open manifolds. Finally, at the discrete 
level, the letter h (with or without geometric meaning) will be used to denote discretiza- 
tion. The expression a < b will be used to mean that there exists C > independent of 
h such that a < Cb for all h. 



2 Variational formulation 

The geometric layout of our problem is as follows: a domain Q G M.'^ {d = 2 or d = 3) 
with polyhedral Lipschitz boundary is subdivided into two subdomains by a Lipschitz 
polyhedral interface E. The subdomains are denoted Qs and ^d (S stands for Stokes and 
D for Darcy). We also denote 

Fs := dns \ E Fd := S^d \ E. 

The normal vector field u on dQ is chosen to point outwards. We also denote by u the 
normal vector on E that points from Qs to Qd. We will add a technical assumption on E 



later on. 

The strong form of the Stokes-Darcy system, as we will consider it here consists of the 
Stokes equations in Qs 

- 2i.div(£(us)) + Vps = fs in ^s, (2.1) 

divus = inf^s, (2.2) 

Us = on Ts, (2.3) 

where s(us) 



^(Vus + (Vug)^), the Darcy equations in fi^, 


K^^ud + VpD = 
divuD = /d 


in fio, 
in fio, 


ud ■ i/ = 


on Td, 






s, conditions 






ud ■ i^ on S 



(2.4) 
(2.5) 

where tt^w := w — (w-i/)i/. The coupling conditions encode mass conservation, balance of 
normal forces and the Beavers- Joseph-Saffman condition. The matrix valued function K 
is assumed to be componentwise in L°°{Qy)), symmetric and uniformly positive definite, 
so that K~^ is componentwise in L°°{Qy)). The coefficient k, G L°°(E) is assumed to be 
bounded from below by a positive constant a.e. on E. Because of the mass conservation 
condition across S, the homogeneous Dirichlet boundary condition for Ug on Fg and the 
incompressibility condition in the Stokes domain, we can easily show that 

/d = (2.6) 

no 

is a necessary condition for existence of solution. 

We set H^(r2g) := H^{QsY and consider the Sobolev spaces 

H^(ng) := {uGHi(fig) : u = on Tg} (2.7) 

HD(div,fiD) := {u G H(div, fio) : u ■ jy = on To}, (2.8) 

endowed with their natural norms || ■ ||Hi(r2s) and || ■ ||H(div,QD) ("^f- II]' [E])- We recall 

I/O 

that the normal trace operator u H- u|s ■ i' is bounded from HD(div, ^d) onto Hqq (S). 
We introduce the space for the velocity field 

X := {u = (ug, ud) G H^(f^g) X HD(div, Qj^) : us ■ v = uj) ■ v on S} C Ho(div, Q), 

(2.9) 
and endow it with the product norm 

,1/2 
|U| 



(l|us|lHi(ns) + l|uD||H(div,nD)j ' V(ug,UD)GX. 



The space for the pressure field is 

Q := Lli^s) X L^(fiD) X M =: L^(fi) x M, 

where 

Ll{0):={peL\0) : {p,l)o = 0}. 

The pressure field is represented as (ps, Pd,S), where ps + S corresponds to the pressure in 
the Stokes domain and normalization has been applied to have zero integral of pressure 
over Qi). If we want to impose zero integral of the pressure field over Q, we can easily 
correct the result by adding a constant in a postprocessing step. The space Q is endowed 
with the corresponding product norm. 

We next define two bounded bilinear forms for the weak formulation of the Darcy- 
Stokes problem. We first consider a : X x X — ?■ M, given by 

a(u, v) := 2z/(£(us), siys))ns + z^(fi;" V^us, 7rtVs)E + (K"^ud, VD)nD- 

We also consider 6 : X x Q — )• M given by 

6(u, (g, 6)) := (div u, g)nsunD + ^("s ■ ^, l)s- (2.10) 

For analytical purposes, this bilinear form can be considered as the sum of two bilinear 
forms 6i : X X Ll{n) ^ M and 62 : X x M ^ M: 

6i(u,g) := (divu,g)nsuCD> b2iu,S) := 5{us ■ u,1)y;. (2.11) 

A variational form of the Darcy-Stokes problem looks for (u, (p, 5)) G X x Q such that 

a(u,v) -6(v, (p,(5)) =(fs,vs)ns Vv G X, 

(2.12) 
6(u, {q, p)) = (/d, qB)nu V(g, p) G Q. 

Terminology. For economy of expression, given a bounded bilinear form c : § x T — ^ M 
(here § and T are Hilbert spaces), we will denote 

kerc:= {s G§ : c(s,t) = Vt G T}. 

(The set is actually the kernel of the operator § 9 s 1— t- c(s, ■) G T*, which will remain 
unnamed.) We will also say that c satisfies the inf-sup condition for the pair {§,¥}, 
whenever there exists /3 > such that 

sup 4^ > P\\th Vt G T. (2.13) 

This property is equivalent to surjectivity of the operator s \-^ c(s, ■ ) and implies the 
existence of a bounded right-inverse of this operator. Actually, given any ^ G T* we can 
find s G S such that 

c(s,t) = e(t) VtGT ||s||s<rieilT*. (2.14) 



The constant (3 in (2.13) and (2.14) can be taken to be the same and the map ^ 1— )■ s is 
linear. 



Proposition 2.1. For the bilinear forms given in (2.10[)-(2.11) it holds that 



ker 61 
kerb 



{uGX 
{uGX 
{uGX 



divuGPo(fis)xPo(fiD)} 
divu G span{|r2t 



1-1- 



- la 



1-1 



^nj} 



divu = 0}, 

where 1^^ and 1q^ are characteristic functions of ^d and Q^ respectively. 

Proof. Note that 

(divUo,go)f7, =0 yqeLlirio) ^^ divUoGPo(fio) og{S,D}. 

Also, for an element u G X 

(divug, 1)^3 = {uu, 1)e = -(divuD, 1)^^- 

The proof of the statement is straightforward with help of these two results. D 

Proposition 2.2. The bilinear form b satisfies the inf-sup condition for the pair {X, Q}. 

Proof. Let 

Xo:=H^(fis)xHo(div,fiD)cX. 

Note first that 61 satisfies the inf-sup condition for the pair {Xo,-L^(fi)}, since this is 
equivalent to separate and well-known inf-sup conditions in fig and ^d- Therefore 61 
satisfies the inf-sup condition for the pair {X, Ll{Q)}. The inf-sup condition of 62 for the 
pair {X, M} is equivalent to the existence of a vector field v G X such that 



(VJ/,1)2>0. 



(2.15) 



In fact, it is possible to construct v G Hj(r2) with this property. The proof of the global 
inf-sup condition for b uses then the characterization given in fiSl Theorem 2]. We need 
to show that X = ker fei -|- ker62- Given v G X we can choose w G Xq such that 



(div w, q)n = (div v, q)n Vg G Ll{n). 



(2.16) 



Note that w = on S and therefore w G ker62 and by definition v — w G kerfei, which 
proves the result. D 



Proposition 2.3. Problem (2.12) is well posed. 



Proof. Sufficient conditions for well-posedness of mixed problems are (see P, Chapter 2]): 

the inf-sup condition of b for {X; i 

kerb 



(shown in Proposition 2.2) and coercivity of a in 



a(u, u) > «||u| 



Vu G ker b. 



(In fact, since a is symmetric and positive definite, these conditions are also necessary.) 
However, the bilinear form a is coercive in the set 



UliUs) X {v G H(div,fiD) : divv = 0}, 



(2.17) 



by Korn's inequality. By Proposition 2.1, this set includes ker 6, which proves the result. 

D 
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3 The discrete problem 

We assume that there are two separate famihes of regular triangulations {T^}h and {T^}h 
of fis and Qb respectively. Each triangulation is composed of triangles/tetrahedra with 
a diameter not greater than h. The triangulations create two inherited partitions of S, 
respectively denoted Eg and S^. Let us consider finite dimensional subspaces of piecewise 
polynomial functions (relatively to the given triangulations) to approximate velocity and 
pressure in the Stokes domain 

U\ns) C H^(fis), L^oi^s) C Llins), L\ns) = Lo"(fis) © Po(fis), 
as well as in the Darcy domain 

U^n^) C H(div, fie), L^{^b) C Llin^), L'^in^,) = L^{^b) © Po(fiD). 
The spaces 

H[;(ns) := H''(fis) n HliQs), Ho"(1^d) := H'^I^d) H Ho(div, Qj,) 

are the corresponding subspaces that are considered when applying the discretization 
method to problems with homogeneous boundary conditions on the entire boundary of 
each subdomain. We also need to consider the spaces (recall (2.7) and ( |2.8[ )) 



H^(fis) := H'^(fis) nH^(f]s), H^(f^D) := H"(f^D) n HD(div,fiD), (3.1) 

as well as the discrete spaces of normal components on S, namely, 

$^ := {u, ■ ^ : u, G H„^(a)} C ^^(S) o G {S, D}. 

In all what follows, we assume that $p contains at least the space of piecewise constant 
functions on E{^, i.e., 

Po(S^):={0/.:S^R : </);,|eGPo(e) Ve G S^} C <l>^. (3.2) 

We denote by R^ the L^(S)-orthogonal projection onto $(^. 

In the forthcoming arguments, we will use the idea of a uniformly bounded right inverse 
of the normal trace H''(r2D) 9 v/^ i-> v/i-i/ G $d, i.e., a linear operator L'' : <I>q — > 11^(^0) 
such that 

||LV"||H(div,nB) < C ll/ll^-;/2(s) V0'^ G <l>^ 

with a constant C > independent of h and 

(lV'') ■iy = (p'' V/ G $1^. (3.3) 

We will refer to such an operator as a stable lifting of the normal trace in H'^(r2D). 

The method we are proposing is a Galerkin discretization of the variational problem 



(2.12) using the spaces 



X'^ := {u, = (ug, u^) G H^(fis) X H^(fiD) : u^ ■ ^ = i?^(u^ u) on S}, 
Q'^ := L^i^s) X 4X^0) X M =: L^(fi) x M, 



that is, we look for (u^, {ph, 6h)) G X^ x Q^ such that 

a{uh, v/,) - b{-Vh, (Ph, Sh) = (fs, v/,)qs Vv,, G X'' 



(3.4) 
b{uh, {qh, ph)) = (/d, qh)nu ^{Qh, ph) G Q''. 



Remark 3.1. Note that^ ^ X unless $3 ^ $d m which case (3.4) becomes a conforming 



Galerkin approximation of (2.12). 



More terminology. The discrete counterpart of the kernel and inf-sup terminology 
introduced in Section [T] is more or less straightforward to define. Let Sh C § and T^ C T 
be sequences of finite dimensional spaces of the Hilbert spaces S and T and consider again 
a bounded bilinear form c : S x T — )■ M. The discrete kernel of c is the set 

kerch := {sh G S^ : c{sh,th) = VU G Th}, 

that is, it is the kernel of the discrete operator Sh 3 Sh ^-^ c{sh, ■ ) G TJ^. We say that 
c satisfies a uniform (discrete) inf-sup condition for the pair {§/i,T/i} when there 
exists /3 > such that 

sup ^^^^>/3\\thh ytneTh yh. (3.5) 

This condition implies that for all ^ G T* we can find a sequence Sh G §h such that 

cish,th) = ^{th) ythETh, \\shh<l3-'m\T* \fh (3.6) 

and the maps ^ ^-^ Sh are linear. 



Let us now discuss the properties of the discrete spaces that ensure stability for (3.4). 
The first set of assumptions is natural for each of the sets of discrete couples. 

Hypothesis 1. (a) The spaces for the discrete Stokes flow are stable when applied to 
the Stokes problem with homogeneous Dirichlet boundary condition: 

(divU/i, g/t)ng ^ o II II w rhfr^ \ m ^\ 

sup \ — , -^>l^s\\qh\\n^ yqheL^ins). (3.7) 

o^UheH{j(ns) l|Uh||Hi(ns) 

(b) The spaces for the discrete Darcy fiow are stable when applied to the Darcy equations 
with homogeneous normal trace: 

(divu/j,gft)no ^ o II II w ^ Th(r> \ ro q\ 

sup - — r^ > PDmhWno Vg/, G Lo(i2D). (3.8) 

Oj^uh&n^^inu) \\'^h\\u{div,Qo) 

(c) divH'^(fiD) cL'^ino). 

Hypothesis 2. (a) There exist Pq, Pi such that for all h 

3v/, G H^(fis) s.t. {vh ■ u, 1)2 > /3o and 1|v;,||hi(Qs) < Pi- (3.9) 



(b) There exists a stable lifting of the normal trace. 
The validity of last set of hypotheses will be discussed in Sections [5] and |6| Although 



at this stage these are just theoretical assumptions of the discrete spaces, condition (3.9) 
can be understood as the possibility of the discrete spaces to create a certain amount of 
flux across E with a bounded velocity field. 



Proposition 3.1. Hypotheses [^ and [^ imply that the bilinear form (2.10) satisfies a 
uniform inf-sup condition for the pair {X'^, Q^}. 

Proof. We will use the schematic method for proving inf-sup conditions in product spaces 
given in p^, Theorem 8]. In order to do it, we consider the decomposition 6 = 6i + 62 



given in (2.11). Let 

xj; := x'^ n Xo = n'^isis) x h^I^d). 



Conditions (3.7) and (3.8) in Hypothesis [T] are equivalent to 61 satisfying a uniform inf-sup 



condition for the pair {Xq,L^(0)}. Therefore, they imply the inf-sup condition of 61 for 

On the other hand, we take Vg as in Hypothesis ^a) and define \^ := L'^i?^(vg ■ u) 
where L^ : $p — )■ 11^(^20) is the stable lifting of the normal trace. Then 



|vS||H(div,n,) < WKiys ■ ^)llH-/^fE^ ^ II^^K ■ ^ 



(2) 



E 



< l|Vs'-^i|E<l 



and Vp ■ fc* = i?^(vg ■ v) which implies that (vg, v^) G X'*. We have then just shown that 
there exist /So, /3i > such that for all h 

3v, = (v|,v^)gX'^ s.t. (v^i/,l)2>/3o and ||v,||x < /?;. (3.10) 



It is easy to see that (3.10) is equivalent to the uniform inf-sup condition of 62 for {X^, M}. 
Let now v/j G X^. By the uniform inf-sup condition of 61 for {Xq, L^{VL)} we can find 
w/j G Xq such that 

(divwft,gft)n = (divv/„g,,)n Vg^, G L^(fi), \\^h\W ^ l|divv;,||n < \\\h\W- 

By definition, w^ G Xq C ker62,?i and v/j — w^ G keihi^h- Therefore we can decompose 
v/i = w/i + (v/i — w/i) G ker 62,/! + ker hi^h, and this decomposition is stable. This is enough 
to prove the uniform inf-sup condition of h. D 

For the sake of measuring errors, we consider the space 

X:=H^(fis)xHD(div,fiD), 

that contains both X and X'^. This space is endowed with the same product norm as 
X. In order to simplify the argument we will use a global bounded bilinear form A : 

(X X Q) X (X X Q) ^ M given by 

^((u, (p, 5), (v, (g, p)) := a(u, v) - 6(v, {p, 6)) + 6(u, (g, p)). 



Proposition 3.2 (Stability and Strang estimate). Hypotheses[^and[^ imply unique solv- 
ability of the discrete equations (3.4) and the error estimate 

||u-Uft,||x+ \\p-ph\\no + \\P + ^ - {Ph + Sh)\\ns 

< inf ||u-Vft||x+ inf \\p - qh\\n + Ch{pD), 

where the consistency error term C/j(pd) =0 if ^^ G $(^ and 

Ch{PD)<h'/'\\pj,-R'},pr,\\j: 
otherwise. 



Proof. We already know from Proposition 3A_ that the inf-sup condition of b holds true 
for the pairs {X'*, Q'^}. Hence, to prove that the operator 

has a uniformly bounded inverse we just have to show that the bilinear form a is coercive 
in the discrete kernel keibh := {u^ G X^ : b{u'^ , {qh, Ph)) = ^{qh-,Ph) G Q'^}. By 
Hypothesis [l[c), if Uh G kerbh, then divu^ G L'^{Q,b) and since 

{diyuh,qh)nu = ^ Wqte Lq{Qb), 

it follows that divu/^ G Po(fiD)- Also, as part of the fact that u^ G kerbh, it follows that 

(divu^, l)f,„ = (u^ ■ ^, l)s = (i?^(ug ■ ly), Ih = i^s ■ ^^ 1)e = 

where we used here the fact that Po(S) C $d- Therefore, if mh G ker6/i, divu/^ = in 
VLy). The coercivity of a on kerfe/^ follows then from the fact that it is coercive on 

H^(fis) X {v G HD(div,fiD) : divv = 0}. 

Let now u := {u,{p,6)) and Uh := {uh,{ph,Sh))- Then, for all Vh = {'Vh,iqh, Ph)) G 
X'^ X Q'^, 

A{u-Uh,Vh) = a(u,v/i) -6(v/„(p, (5)) - (fs,v/i)ns (3.11) 

and 

a(u,v) -6(v,(p, (5)) - (fs,vs)ns = -(vs -J^ - vd ■j^,Pd)e Vv G X. (3.12) 

In this last formula, we have used that if (u, {p,S)) G X x Q is the solution of the equations 
( |2l2| , thenpD G //^(^d)- Therefore, by ( [37II| and ( |3l2| , 

A{u-Uh,Vh) = -{{^r|-^r'^)■ly,PJ,)s. (3.13) 

The discrete inf-sup condition satisfied by the global bilinear form A and ( 3.13[ ) show 
that for all w^ G X'' x 



h ^ inh 



\u-Uh\\xxQ < \\u-Wh\\xxQ+\\Uh-Wh\\xxQ 

A{uh - Wh,Vh) 



< \\u-Wh\\xxQ+ sup 

^ II- - II , A{u-Uh,Vh) 

< \\u-Wh\\xxQ+ sup p— ^ 

OT^S-^eX'^xQ'' ll'i^/illXxQ 

I((vs'-v^)-^,Pd)^ 



< 



\u-Wh\\xxQ+ sup 



o^irhex''xQh |p/i||xx(! 
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It is clear that, if $g C $^5 ^^e last term of the previous inequality vanishes identically. 
In the general case we have to estimate the consistency error 

a.(po):^ sup K(v§-vfe)-.Po)d 



To this end we introduce the L^(S)-projection IIq : i^^(S) — )■ Po(Sq) onto the space of 
piecewise constant functions and denote IIq : L^(S)'^ — ;■ Po(Sq)'^ its vectorial counterpart. 
It is straightforward that 

((Vs-VD)-i^,PD)E = ((Vs-v{^)-jy,pD-i?DPD)s = (Vs-J^,PD-i?DPD)s 

= (Vs^■^-<(v^^),pD-i?^J9D)E 

where in the last inequality we have applied a well known approximation estimate for 
piecewise constant functions. Using the trace theorem in H^(r2g) we deduce that the 
consistency error may be bounded by 

ChiPD) < h'/^pj, - R'hPoh 

and the result follows. D 



4 Approximation properties of 



'h 



In principle, the restriction of equal normal component on S can reduce the size of the 
separate spaces (H^(f2s) and H^^Qb)) and limit the approximation properties, unless 
appropriate matching on the boundary allows for a rich enough discrete space. Note that 



as proven in Proposition 3.2), this is not a matter of stability, but of approximation. 



Proposition 4.1. Hypothesis\^b) implies that for a// u G X 



inf ||u-v,,||x< inf ||ug - Us||hi(Qs) + i^f ||ud - UD||H(div.QD) 
VhSX'' ugeH|{Qs) u£eH£(nD) 

+ A(/i)||us-i/-i?^(us-iv)||E, (4.1) 

with \{h) =0 z/<l>g C $0 ^'^^ K^) ~ ^^^^ otherwise. 
Proof. Let 

n^Hi(fis)-^H|(fis) and H^ : Holdiv, ^d) ^ H^I^d) 



be the orthogonal projections onto the two discrete spaces (3.1 ). Given u = (us, Uq) G X, 
we consider 

u, = (u^,ui^) := (n^us,n^uD - L'^(n^uD ■ ly - i?^(n^us ■ i^))). 
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Since u^ ■ u = i?^(ngUs ■ u) = R^i^s ■ t^), it follows that u/j G X^. On the other hand, 
the uniform boundedness of L^ yields 

l|u - Uftllx < ||us - ngUsllHi(f^s) + II^D - nDUDl|H(div,QD) + 

iKuD ■ ly - /2^(n^us ■ ^)IIh-/^(e) (4-2) 
and the triangle inequality together with the fact that ug ■ v = ud ■ iy gives 

llnJ^UD ■ ly - RUU^ns ■ ^) 11^-1/2(2) < II (ud - n^uo) • ^^11^-1/2(2) + 

||us-^-i?^(n^us-^)||^-i/2(^). (4.3) 



Consequently, if $g C $u, -R^ is redundant in the last estimate. Consequently, (4.1) 
is satisfied in this case with A(/i) = thanks to the boundedness of the normal trace 
operator. 

In the general case, we will need the estimate 

U-RU\\H-,^f^^^)^h'^"Uh VeGL2(S), (4.4) 



obtained from a duality argument by taking into account (3.2). We estimate the second 



term of the right hand side of (4.3) by the triangle inequality 



|us • ^ - Rni^s^s ■ ^)IIh,-/''2(s) ^ ll"s ■ ^ - Rni^s ■ J^)\\h-,'/\j:) + 

||US • I' - HgUg • i^||^-l/2(2) + ll(id - ^d)(uS ■ 1^ - ngUg ■ '^)||^-i/2(s) 



and by using (4.4) twice 

||us • ly - i?^(n^us ■ ^)ll^j/2(s) ^ ^'^' ll"s ■ ^ - RU^s ■ ^)||e + ||us - n^ug||s. 
The result follows by applying the trace theorem in H^(r2s) and combining the resulting 



estimate with (4.2). D 



We are now in a position to establish our main result. 



Theorem 4.2. Hypotheses^ and^imply unique solvability of the discrete equations (3.4) 
and the error estimate 



\u-Uh\\x+ \\p-Ph\\no + lb + 5- {Ph + 5h)\\ns 



< 



inf ||ug-u^||Hi(ns) + inf ||ud - u|^||H(div,nD) + inf \\p - QhWa 

+ A(/i) (Ibo - R'hPBh + ll"s ■ ^ - RU^s ■ ^)||s) 
where X{h) = z/$g C $d ^'^^ H^) ~ ^^^"^ otherwise. 
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5 The lifting of the normal trace 



Theorem 4.2 shows that, given a set of two mixed finite elements that are separately con- 
vergent for the Darcy and the Stokes problems, the only requirement for the convergence 
of scheme ( |3.4[ ) is Hypothesis 12} In this section, we discuss the conditions under which 
Hypothesis |2[b) is satisfied for the most common finite element subspaces of H(div, ^d). 
We begin by proving the existence of a stable lifting of the normal trace for Raviart- 
Thomas and Brezzi-Douglas-Marini spaces in two dimensions. The only restriction on the 
grid is shape-regularity. For the sake of clarity, we define all needed elements here. The 
geometric elements are a polygonal domain fi C M^ with boundary F, a shape-regular 
family of triangulations Th , and the partitions F^ of F that are inherited from Th- The 
space of vector discrete vector fields is the following: 

H'^ := {v;, e H(div,fi) : v|t G P(T) VT G TJ, 

where 

^ ' •" \ Pfc(T)2, if A; > 1. 

Note that H'' is the velocity space for the Raviart-Thomas pair when k = and for the 
BDM pair for k > 1. The latter is a strict subspace of the RT space of the same order. 
The space of normal traces is 

^'^-{aiF^M : aieGPfc(e) Ve G F^. 

It is clear that u^ ■ i/ E ^h for all Uh G H^. 

Theorem 5.1. There exists L'* : ^'^ — )• H'^ such that 

(L'^a)-^ = a and ||L'^al|H(div,f^) < ||all/f-i/2(r) V^ G <l'^ (5.1) 



In addition, 



divL'^a = 7^ /a- (5.2) 



Proof. We start by decomposing 

a = ^Jjh+ {^h- ^Jjh) = c, + eke Po(F) © ($, n h,'/\t)) , (5.3) 

where H^'^\T) := {^ G H-^^T) : {^, l)r = 0}. 

Lifting of a constant function. Consider the solution of the Neumann problem 



A?i = \n\'^ in n, d^u = \T\-\ u = 

Jn 

and let w := Vm. From well known regularity results [15] it follows that w G H^/^+^(i7) 
for some e = s{Q) > 0. Let then Wh be the lowest order Raviart-Thomas projection of 
w, i.e., 

w;,GH(div,fi), w;,GPo(T)2©span{x} WT e %, / w^ ■ zv^ = / w • i/^ Ve G ^fe, 
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where £h is the set of edges of the triangulation and Ue is the unit normal vector on e 
(for any given orientation). This projection is well defined because of the regularity of w. 
It is then clear (using a direct argument or well known properties of the Raviart-Thomas 
element) that 

w,,-z/=|r|-\ divw/,= |f]|"^ (5.4) 

A more delicate argument (see for instance pj) shows that 



|w/i||H(div,Q) ^ ||w||H(div,n) + ||w||Hi/2+e(n) — : C'sn 



(5.5) 



Lifting by arc-length integration. The gist of the lifting process follows from the application 
of a projection on a space of continuous finite elements after integration of ^^ along the 
length of r. Consider the spaces 



V'' := {UH E C(fi) : Uh\t G Pfc+i(r) VT G %}, 



•q/^ 



^V^ = {^Uh : Uh e V''} 



where 7 stands here for the trace operator on F. In [21] an operator Sh '■ H^{Q) — t- V^ is 
constructed with the following properties: 



Slu 



ShU, \\Shu\\HHn) < \\u\\m(n) Vm g H^i^), 



and 



7M = 



'jShU = 0. 



(5.6) 



(5.7) 



Note that (5.6) and (5.7) imply that if 7U G \I' , then ^yShU = ^u. (This property can be 



verified directly for the particular construction of Scott and Zhang [21], or proved directly 
from properties (5.6) and (5.7).) Let then D^^ : H^ (F) — )■ H^/'^{T) be an inverse 
of the arc-length differentiation operator and let L : H^^'^iV) -^ H^iVt) be a bounded 
right-inverse of the trace operator. Given ^° G $'' fl /Jq (r), we can define 



V, := y^S^LD-^il V^ := {dy, -d^). 

Note that -fLD'^^l = D-i^° G *^ and therefore -fShLD-^^^ = D'^^l Hence 



1-ltON/ 



V, ■ z. = i^S^LD-'eh)' = {D-X 






{5.t 



It is also clear that \h is piecewise Pfc and that its divergence vanishes. Therefore v/j G H'^. 
Finally, 



|V/i||H(div,n) 



|vft||n < \\ShLD ^^°||//i(n) < \\LD ^^°||/^i(n) 



< \\D''eH\\mr^ir) < UhWH-^/Hry (5.9) 
Conclusion. The full lifting operator is then given by the expression 



using the discrete lifting of a constant given above (5.4)-(5.5). By (5.4) and (5.8), it 
follows that 



(L'^a) ■ ly = ^Jjh+ [^h- ^Jjh) = Ch 
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Also, by (5.5) and (5.9), it follows that 



Lf^c II <^ 

, ?/i||H(div,n) ^ 



a 



IK/i||//-i/2(r) < IK?»||/f-i/2(r)- 



Finally, by (5.4) 



divL'^e/. 



^h ) div Wh 



\n\ 



which finishes the proof. 



D 



Remark 5.1. In the three dimensional case a stable lifting of the normal trace can be 
constructed for RT and BDM elements using an additional hypothesis on Th, namely, that 
Th is quasiuniform on a neighborhood ofV . A proof for RT of any order in two dimensions 
is given in fll^ : the proof holds for three dimensions and for BDM elements. 



6 The discrete permeability condition 

In this section we discuss if Hypothesis [2|a) is available in practical situations. This 
condition can be understood as the possibility of having fiow across S without the need of 
an increasingly large velocity field, that is, the discrete boundary created by the choice of 
spaces has to be permeable enough. Let us first show that Hypothesis |2](a) can be easily 
deduced from Hypothesis [2| 



Proposition 6.1. Assume that Hypothesis^b) is satisfied and that 

uj{h) := inf ||us - Ug||Hi(ns) + inf ||ud - UD||H(div,f7D) 

converges to zero when the parameter h tends to zero. Then, there exists Hq > such that 
Hypothesis IWa) is satisfied for all h < h^ 

Proof. Let u G HQ(fi) be such that 

||u||x < 1 and (u ■ i/, l)s > 1. 

With the notations of Section H] we define 

u, ^ (u^, O := (n^us, n^uD - L'^(n^uD ■ ^ - i?^(n^us ■ ^))), 



where Ud := uJ^d and us := uj^g. We know from Proposition 4.1 that 



l|u - ^h\\x ^ ^(h) + \{h) ||us ■ u - R^ius ■ J^)||s, 
with \{h) < h^/'^. It follows that 

||u/i||x < \\uh\\x + ||u - Uh\\x < 2 

and 

{uh ■ V, l)s = {uu, 1)e - ((u - Uh) • J^, l)s > 1 - ||u - UhWx > 1/2 

if h is sufficiently small. 



D 
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Next, we will see that Hypothesis |2Va) is actually quite mild, it can be shown to hold on 
(independently from Hypothesis [2Fb)) with few restrictions. In some of the forthcoming 
arguments and examples we will use the space of continuous finite elements: 

PTiri") := {v G CiUsY : v|t E Pfc(T) VT E Tg"}. 
Proposition 6.2. In the two dimensional case, assume that: 

(a) either H is a polygon with two or more edges, 

(b) orHisa line segment and there exists a fixed node p belonging to all partitions Eg 
such that the distance of p to dVt remains hounded below. 

Additionally, assume that 

pcont ^j-h-^ C H'*(fis). (6.1) 

Then Hypothesis^a) is satisfied. 

Proof. We start with the geometric condition (a). We choose two adjacent edges. Si and 
S2, of S and let p := Si fl S2. We next construct the function p : S — )• M^ given by 

where Ui is the normal vector on Sj and ■?/' : S — )■ M is a continuous piecewise linear 
function (piecewise with respect to the natural partition of S in edges and linear with 
respect to the arc parameterization) such that V'(p) = 1 ^ind if) vanishes on all other 
vertices of S. Note first that 

1 /" 1 

pu = - (1 + ^/1-^/2)^= t(1 + J^i-J^2)|SiUS2| =: Co >0 (6.2) 

S ^ JS1US2 4 

with Co independent of any discrete quantity. We next use the interpolation operator for 
non-smooth functions of Scott and Zhang ^2lj (see also |1] and [^) to construct v/^ E 
Pi°°*(7^'^) such that v/i = p on S and v/^ = on Fg. Because the Scott-Zhang lifting 
operator is stable, it follows that 



■'oo 



where Co independent of any discrete quantities. Since (6.1) is satisfied, this proves the 
hypothesis. 

In case (b), we can choose Si and S2 to be the segments joining p to the boundary 
of VL and proceed similarly. In this case, the quantity Cq in ( 6.3[ ) depends on p: the 



hypothesis on the distance of p to dVt is equivalent to having min{|Si|, IS2I} > Ci > 0, 
which can be used to give an upper bound of ||p||„i/2,„,. D 



Remark 6.1. The idea of Proposition 6.2 can he extended to cover more cases: 

(a) If¥^^^^{T^) C H'^(fis) (o-s usual, d is the dimension of the physical space), a sim- 
plified proof can he carried out hy choosing a face of S, huilding a huhhle function p 
in the normal direction and extending it to f2g using the Scott-Zhang interpolation 
method. If S does not contain any triangular face, we have to additionally assume 
there exists a fixed triangle, that can always he ohtained as the union of elements 
(triangles) from Sg. 
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(b) In the three dimensional case, with lower order polynomials, it is easy to find geo- 



metric configurations of E that match the requirements of Proposition 6^, so that 
we can construct a discrete bubble function on an averaged normal direction on part 
of the interface. 

7 Examples 

In this section we provide several examples of pairs of stable elements for the Stokes-Darcy 
problem. For careful description of the associated mixed elements on the Darcy domain 
and for stable finite elements for the Stokes problem, the reader is referred to [5j and 
[9]. Original sources for these elements (many of them discovered in several steps and 
renamed several times) can be found in these already classical references. Spaces in the 
Stokes and Darcy domain will be chosen so that Hypothesis [T] is satisfied. 

7.1 The conforming case 

We next list several examples of choices where the hypotheses above are met. In all the 
examples below, the space for the Darcy domain will be the Brezzi-Douglas-Marini (BDM) 
space, sometimes referred to as the Brezzi-Douglas-Duran-Fortin in the three dimensional 
case. For a triangulation 7^ of ^d, we consider the spaces for k > 1: 

H\Q^) := {u, GH(div,fiD) : u,|rGPfc(T)'^ ^TeT^}, 
L'^(nD) := {Ph-^D^^: PhlreFk-iiT) ^T eT^}. 

We will refer to this choice with the generic name BDM(/i;). If the inherited triangulation 
of E is denoted E^, the space $q consists of piecewise P^ functions. For conditions on 
when there is a stable discrete lifting, the reader is referred to Section [5j 

In the Stokes domain, we consider another triangulation Tg^, producing a partition 
Eg of the interface. We assume that the Darcy partition E^ is either equal to or a 
refinement of Eg. Since all discrete spaces for the velocity variable in the Stokes domain 
contain polynomial linear functions, we are going to assume that the geometry allows for 
Hypothesis [2|a) to be satisfied. 



A list of possible choices is given in Table 7.1 



7.2 The nonconforming case 

In addition to the BDM element, we now consider the Raviart-Thomas element (also 
called Raviart-Thomas- Nedelec in the three dimensional case). The spaces for the RT(A;) 
pair with k > are 

ll\n^) := {u;, e H(div, Qd) : u;,|t G niT)" © Pfc(r){x} VT G T^}, 

' v/ ' 

=RTfc(T) 

L\n,,) := {ph-.Qb^^: PhWeMT) \fTeT^}. 

For this kind of coupling the Stokes and Darcy triangulations can be taken to be com- 
pletely independent (although this might seriously complicate implementation in the three 
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Stokes 


Velocity 


Press. 


Darcy 


Vel. 


Press. 


Order 


MINI 


Pi+bubbles 


pcont 


BDM(l) 


Pi 


Po 


h 


Taylor-Hood, k>2 


Pfc 


"pcont 


BDM(fc) 


Pfc 


fk-i 


h^ 


Conf Crouzeix-Raviart 


Ps+bubbles 


Pi 


BDM(2) 


P2 


Pi 


h' 


Bernardi-Raugel 


Pi+face bubbles 


Fo 


BDM(l) 


Pi 


Po 


h 



Table 1: Coupling of Stokes elements with BDM elements. The superscript '^°'^* refers 
to the demand of continuity for the discrete pressure space. The bubbles are used for 
velocities in the MINI and conformal CR elements: an internal Fd+i{T) bubble is added 
to the velocity space on each element. For the BR element, face bubbles are included on 
all internal faces, but no bubbles are added on faces lying on E. When these bubbles (no 
needed for stability) are added, the method stops being a particular case of this class. 



dimensional case). Condition (3.2) is satisfied by the IiT{k) spaces A; > and the BDM(/i;) 
spaces k > 1. For conditions concerning the lifting of the normal trace, see Section |5} 
All Stokes spaces contain piecewise linear functions, which means that except in trivial 
geometric configurations Hypothesis ^a) will be satisfied. 



Stokes 


Velocity 


Press. 


Darcy 


Vel. 


Press. 


Order 


MINI 


Pi+bubbles 


pcont 


RT(0) 


RTo 


Po 


h 


Taylor-Hood, k>2 


Pfc 


■pcont 
^fc-1 


RT(A; - 1) 


RTfe-i 


Pfc-i 


h'' 


Bernardi-Raugel 


Pi -|- face bubbles 


Po 


RT(0) 


RTo 


Po 


h 


Ps-iso-Pi 


P^(rV2) 


pcont 


BDM(l) 


Pi 


Po 


h 



Table 2: Coupling of Stokes elements with BDM and RT elements and their order of 
convergence. The superscript ™°* refers to the demand of continuity for the discrete 
pressure space. The bubbles are used for velocities in the MINI element. For the BR 
element, face bubbles are only included on the internal faces. Adding them to faces on E 
does not change the convergence order. In that case BR can be coupled with BDM(l) as 
well. 



8 Numerical results 



In order to confirm the good performance of our scheme (3.4), we present in this section the 
combination of several stable Stokes elements with the Raviart- Thomas element and the 



Brezzi-Douglas-Marini element, as shown in Tables 7.1 and 7.2 We begin by introducing 
some notations. The variable N stands for the total number of degrees of freedom defining 
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the finite element subspaces X'^ and Q'^, and the individual errors are denoted by: 



e(UDJ :— ||ud — Uj-,||H(div,CD)5 



eiusj := ||us - Ug||Hi(Cs); 



and 



e(pD) := Wpd - Puhu, e(ps) := lbs - Pshs, 

where u'^ := Uh\no, Ug := u^^l^^,, p'^ := ph\no and pg := ph\ns + ^h with (u/^, {ph,Sh)) G 



j^h ^ Qh ]3gi]2g the solution of (3.4). We also let r(uD), r(us), r{pD) and r{ps) be the 

log(e(us)/e'(us)) 
and 



experimental rates of convergence given by 

log(e(uD)/e'(uD)) 



log(e(pD)/e'(pD)) . ^ 

'"(Pd) := i — 77:7777 , r{ps) 



\og{h/h') 
log(e(ps)/e'(ps)) 



log(/i//i') ' '"^^^ ■ log(/i//i') 

where h and /i' are two consecutive mesh sizes with errors e and e'. 

Figure 1: MINI-BDM(l) coupling (left) and Bernardi-Raugel-RT(O) coupling (right). 





Table 3: Convergence rates. 





MINI-BDM(l) 


MINI-RT(O) 


hs 


/id 


''(us) 


''(ud) 


rips) 


r{pB) 


''(us) 


'■(ud) 


r{ps) 


ripu) 


1/6 


1/12 


- 


- 


- 


- 


- 


- 


- 


- 


1/10 


1/20 


1.091 


0.983 


1.990 


1.017 


1.091 


0.979 


1.990 


0.988 


1/14 


1/28 


1.063 


0.993 


1.955 


1.010 


1.063 


0.991 


1.955 


0.995 


1/18 


1/36 


1.047 


0.996 


1.926 


1.006 


1.047 


0.995 


1.926 


0.997 



We now describe the data of the example. We consider the domains Qb := (0, 1)^ x 
(0,0.5)^ and Qs '■= (0, 1)^ x (0.5, 1), and take z/ = 1, k = 1 and K = I, the identity of 
]^3x3_ Non-homogeneous transmission conditions are considered in order to have an exact 
solution given by: 

puix) := si(l - xi) sin(27ra;i) 2:2(1 - X2) sin(27ra;2) X3 sin(27rx3) - poo , 



19 



Table 4: Convergence rates. 





MINI-BDM(l) 


MINI-RT(O) 


hs 


hu 


''(us) 


r-(uD) 


rips) 


"r{VT>) 


r-(us) 


r(uD) 


r{vs) 


r{VT,) 


1/12 


1/6 


- 


- 


- 


- 


- 


- 


- 


- 


1/20 


1/10 


1.063 


0.933 


1.907 


0.996 


1.061 


0.921 


1.907 


0.951 


1/28 


1/14 


1.037 


0.971 


1.856 


1.021 


1.037 


0.966 


1.856 


0.980 


1/36 


1/18 


1.027 


0.984 


1.824 


1.018 


1.026 


0.981 


1.824 


0.989 



in the porous media and by 

-2xi(l-xi)(l-2x2)(l-2x3) 
us(a3) := xi(\ - xx)x'2,(\ - X2)x-i{\ - 0:3) | x^il - x-2){\ - 2xi)(l - 2x3) 

X3(l-a;3)(l-2xi)(l-2x2) 



and 

in the Stokes domain. 



ps(a;) := exp(xi + X2 + X3) 



Table 5: Convergence rates. 





Bcrnardi-Raugel- 


-BDM(l) 




Bcrnardi-Raugel 


-RT(0) 




hs 


/iD 


r(us) r(uD) 


r{vs) 


"r{pT>) 


r(us) 


r(uD) 


r{vs) 


AVT,) 


1/6 


1/12 


- - 


- 


- 


- 


- 


- 


- 


1/10 


1/20 


1.537 0.983 


1.034 


1.020 


1.537 


0.979 


1.038 


0.991 


1/14 


1/28 


1.574 0.993 


1.021 


1.011 


1.574 


0.991 


1.021 


0.996 


1/18 


1/36 


1.574 0.996 


1.014 


1.007 


1.574 


0.995 


1.014 


0.998 



Table 6: Convergence rates. 





Bcrnardi-Raugel- 


-BDM(l) 




Bcrnai 


di-Raugcl 


-RT(0) 




H 


hu 


r(us) r(uD) 


r{ps) 


r{PD) 


^(us) 


''(ud) 


r{ps) 


r{PD) 


\lvi 


1/6 


- - 


- 


- 


- 


- 


- 


- 


1/20 


1/10 


1.573 0.933 


1.014 


0.996 


1.573 


0.921 


1.014 


0.951 


1/28 


1/14 


1.543 0.971 


1.008 


1.021 


1.543 


0.966 


1.008 


0.980 


1/36 


1/18 


1.503 0.984 


1.005 


1.018 


1.503 


0.981 


1.005 


0.989 



The numerical results were obtained using a MATLAB code. In Figure[T]we summarize 
the convergence history of the Galerkin scheme (3.4) for a sequence of uniform meshes of 
the computational domain Q := (0, 1)^ by means of tetrahedra. We select the conforming 
example consisting in the MINI-BDM(l) coupling and the nonconforming example given 
by the Bernardi-Raugel-RT(O) coupling. In each case we display the individual errors 
versus the degrees of freedom A^. We observe that, as expected, the convergence is linear 
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with respect to the discretization parameter h in all the unknowns unless for the fluid 
pressure in the MINI-BDM(l) example where a quadratic convergence is attained. We 
notice that the Bernardi-Raugel-RT(O) case delivers a convergence in the fluid velocity 
that is slightly faster than 0{h). 

We also provide numerical results for triangulations with hanging nodes on the trans- 
mission interface. We consider, uniform triangulations of the subdomains Qn := (0, 1/2)^ 
and Qs '■= (1/2, 1)^ with a mesh size in one of the subdomains equal to half the mesh size 
in the other one. The expected rates of convergence are attained in all the (conforming 
and non-conforming) cases considered through Tables |3| |4| [s] and |6} 

Summarizing, the numerical results presented here constitute enough support to our 
theory for the strong mixed finite element coupling of Darcy-Stokes flow problem. In a 
forthcoming work we will discuss an efficient iterative method to solve the linear system 
of equations arising from our discretization method. 
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